数学模型

Wang Haihua

🍈 🍉🍊 🍋 🍌

利用均匀分布的随机数可以产生具有任意分布的随机变量的样本，从而可以对随机变量的取值情况进行模拟。

离散随机变量的模拟

设随机变量 $X$ 的分布律为 $P\left\{X=x_{i}\right\}=p_{i}, i=1,2, \cdots$, 令 $p^{(0)}=\mathbf{0}, \quad p^{(i)}=\sum_{j=1}^{i} p_{j}, i=1,2, \cdots$, 将 $p^{(i)}$ 作为分点, 将区间 $(0,1)$ 分为一系列小区间 $\left(p^{(i-1)}, p^{(i)}\right)(i=1,2, \cdots)$ 。对于 均匀分布的随机变量 $R \sim U(0,1)$, 则有 $$ P\left\{p^{(i-1)}<R \leq p^{(i)}\right\}=p^{(i)}-p^{(i-1)}=p_{i}, \quad i=1,2, \cdots, $$ 由此可知, 事件 $\left\{p^{(i-1)}<R \leq p^{(i)}\right\}$ 和事件 $\left\{X=\boldsymbol{x}_{i}\right\}$ 有相同的发生概率。因此我 们可以用随机变量 $R$ 落在小区间内的情况来模拟离散的随机变量 $X$ 的取值 情况。

问题1

已知在一次随机试验中, 事件 $A, B, C$ 发生的概率分别为 $0.2$, $0.3,0.5$ 。模拟 100000 次随机试验, 计算事件 $A, B, C$ 发生的频率。 在一次随机试验中, 事件发生的概率分布见表1。

表 1 事件发生的概率分布
事件 A B C
概率 0.2 0.3 0.5
累积概率 0.2 0.5 1

分析

我们用产生 $[0,1]$ 区间上均匀分布的随机数, 来模拟事件 $A, B, C$ 的发生。由上述数据和几何概率的知识, 可以认为如果产生的随 机数在区间 $[0,0.2]$ 上, 事件 $A$ 发生了; 产生的随机数在区间 $(0.2,0.5]$ 上, 事件 $B$ 发生了; 产生的随机数在区间 $(0.5,1]$ 上, 事件 $C$ 发生了。产生 100000 个 $[0,1]$ 区间上均匀分布的随机数, 统计随机数落在相应区间上的次数, 就是在 这 100000 次模拟中事件 $A, B, C$ 发生的次数,再除以总的试验次数 100000 , 即得到事件 $A, B, C$ 发生的频率。

代码

from numpy.random import rand
import numpy as np
n=100000; a=rand(n); n1=np.sum(a<=0.2)
n2=np.sum((a>0.2) & (a<=0.5)); n3=np.sum(a>0.5)
f=np.array([n1,n2,n3])/n; print(f)

问题2

事件 $A_{i}(i=1,2, \cdots, 10)$ 发生的概率如下表所示, 求在 10000 次模拟中, 事件 $A_{i}(i=1,2, \cdots, 10)$ 发生的频数。

事件 A_(1) A_(2) A_(3) A_(4) A_(5) A_(6) A_(7) A_(8) A_(9) A_(10)
概率 0.2 0.05 0.01 0.06 0.08 0.1 0.3 0.05 0.03 0.12

代码

from numpy.random import rand
import numpy as np
n=10000; a=rand(n);
p=np.array([0.2,0.05,0.01,0.06,0.08,0.1,0.3,0.05,0.03,0.12])
cp=np.cumsum(p); c=[]; c.append(np.sum(a<=cp[0]))
for i in range(1,len(p)):
    c.append(np.sum((a>cp[i-1]) & (a<=cp[i])))
print(c)

连续随机变量的模拟

利用 $[0,1] 区$ 间上的均匀分布随机数可以产生具有给定分布的随机变量数列。

我们知道, 若随机变量 $\x_i$ 的概率密度函数和分布函数分布为 $f(x), F(x)$ 则随机变量 $\eta=F(\x_i)$ 的分布就是区间 $[0,1]$ 上的均匀分布。因此, 若 $R_{i}$ 是 $[0,1]$ 中均匀分布的随机数, 那么方程 $$ \int_{-\infty}^{x_{i}} f(x) d x=R_{i} $$ 的解 $x_{i}$ 就是所求的具有概率密度函数为 $f(x)$ 的随机抽样。这可简单解释如下。

若某个连续型随机变量 $\xi$ 的分布函数为 $$ F(x)=\int_{-\infty}^{x} f(t) d t, $$ 不失一般性, 设 $F(x)$ 是严格单调增函数, 存在反函数 $x=F^{-1}(y)$, 下面我们 证明随机变量 $\eta=F(\xi)$ 服从 $[0,1]$ 上的均匀分布, 记 $\eta$ 的分布函数为 $G(y)$, 由 于 $F(x)$ 是分布函数, 它的取值在 $[0,1]$ 上, 从而当 $0<y<1$ 时 $$ G(y)=P\{\eta \leq y\}=P\{F(\xi) \leq y\}=P\left\{\xi \leq F^{-1}(y)\right\}=F\left(F^{-1}(y)\right)=y, $$

因而 $\eta$ 的分布函数为 $$ G(y)= \begin{cases}0, & y \leq 0 \\ y, & 0<y<1 \\ 1, & y \geq 1\end{cases} $$ $\eta$ 服从 $[\mathbf{0}, 1]$ 上的均匀分布。 $R$ 为 $[0,1]$ 区间上均匀分布的随机变量, 则随机变量 $\xi=F^{-1}(R)$ 的分布函数 为 $F(x)$, 概率密度函数为 $f(x)$, 这里 $F^{-1}(x)$ 是 $F(x)$ 的反函数。

所以, 只要分布函数 $F(x)$ 的反函数 $F^{-1}(x)$ 存在, 由 $[0,1]$ 区间上均匀分布 的随机数 $R_{t}$, 求 $x_{t}=F^{-1}\left(R_{t}\right)$, 即解方程 $$ F\left(x_{t}\right)=R_{t}, $$ 就可得到分布函数为 $F(x)$ 的随机抽样 $x_{t}$ 。

问题3

求具有指数分布 $$ f(x)= \begin{cases}\lambda e^{-\lambda x}, & x>0 \\ 0, & x \leq 0\end{cases} $$ 的随机抽样。

建模

设 $R_{i}$ 是 $[0,1]$ 区间上均匀分布的随机数, 得 $$ R_{i}=\int_{-\infty}^{x_{i}} f(x) d x=\int_{0}^{x_{i}} \lambda e^{-\lambda x} d x=1-e^{-\lambda x_{i}} . $$ 所以 $$ x_{i}=-\frac{1}{\lambda} \ln \left(1-R_{i}\right) $$ 就是所求的随机抽样。

由于 $1-R_{i}$ 也服从均匀分布, 所以上式又可简化为 $$ x_{i}=-\frac{1}{\lambda} \ln R_{i} . $$

Python 中 NumPy 库函数可以产生常用分布的随机数, 实际上不需要 我们去生成各种分布的随机数。